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ABSTRACT 

In a cold-dark-niatter universe, cosmological structure formation proceeds in rough 
analogy to origami folding. Dark matter occupies a three-dimensional 'sheet,' non- 
intersecting in six-dimensional velocity-position phase space. At early times, the sheet 
was flat like an origami sheet, i.e. velocities were essentially zero. The present paper 
further illustrates this analogy, and identifies a result of origami mathematics that 
applies to cosmology. We define caustics in the initial conditions (Lagrangian space) 
as surfaces in this sheet, along which the sheet has folded. The regions outlined by 
these caustics, which we call streams, may be colored according to the two possible 
orientations of initial basis vectors, a two-coloring such that adjacent streams are not 
colored the same. While this may not have clear observational consequences, it is a 
severe restriction on connectivity, since there are no bounds on the number of colors 
required to color a general arrangement of three-dimensional regions. Then, measuring 
the relevant quantities from A'^-body simulations, we explore how well outer caustics in 
Lagrangian space correspond to a Zel'dovich prediction, as well as to a measurement 
from the recent ORIGAMI algorithm. 
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1 INTRODUCTION 

On moderately large scales (not so large that it is entirely 
described by linear theory), the matter and galaxies in the 



Universe trace what is known as a cosmic web (e.g., Bond 



|et al. 19961. Regions well under the mean density develop 



into voids. At higher density, the matter is arranged into 
planar structures known as walls or pancakes. Still higher- 
density matter has a filamentary morphology. At the highest 
densities, it is assembled into pointlike clusters, or haloes. 

There are several ways to characterize these struc- 
tures. Observationally, one way of identifying and defining 
them is by looking for depressions, ridges, and peaks in 
the (over)density field 5. For example, one can measure the 
eigenvalues of the Hessian d^5 / dxidxj; this gives local den- 
sity maxima/ridges in one, two and three dimensions (e.g., 
Hahn et al.||2007| |Arag6n-Calvo eFaL][2007l [Sousbie et aL] 



2008 



Sousbie||2011| ). Another approach is global, defining 



voids to tessellate space. For instance, they can be defined 
as density depressions outlined by a watershed transform 



(Platen et al. 2007 Neyrinck 2008 1. In this framework, walls. 



filaments, and haloes are defined according to where voids 
meet each other, and the dimensionality of borders separat- 
ing them ( Aragon-Calvo et al. 2010|. In cosmology, these 
two definitions turn out to be rather similar, but this only 



is true in detail if no locally defined walls and filaments end 
within voids. 

Another way to understand the structures is as folds in a 
three-dimensional manifold ('sheet'). This is most suited to a 
Lagrangian viewpoint, in which mass parcels are tracked. A 
complementary viewpoint is called Eulerian, in which densi- 
ties and velocities at fixed positions are tracked. First-order 
Lagrangian perturbation theory, known as the | Zel'dovich] 
( |1970[ ) approximation, already is quite useful to understand 
the basics of cosmological structure formation. 

In this Lagrangian framework, particles in an TV-body 
simulation can be thought of not simply as blobs of mass, but 
as vertices on an initially cubic grid that gravity deforms, 
and eventually causes to self-cross in three-dimensional posi- 
tion space. In six-dimensional velocity-position phase space, 
though, the manifold does not cross itself, a property which 
prompted our recent analogy to the paper-folding that oc- 
curs in origami (Falck et al. 2012 1. In that paper, we use 



an algorithm called ORIGAMI (Order-Reversing Gravity, Ap- 
prehended Mangling Indices) to identify structures in an 
A-body simulation according to the number of axes along 
which the initial-conditions (Lagrangian) lattice has crossed 
itself. A couple of other recent papers have also explored 
the power of comparing initial to final conditions explicitly 



(Shandarin et al. 2011 Abel et al. 20111, for the purposes 
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of finding shell-crossings, and measuring the density within 
the phase-space sheet. 

Another major application of a Lagrangian framework 
is the study of caustics, largely the subject of present work. 
A 'caustic' is the edge of a fold of the sheet in six-dimensional 
phase space. A strong observational motivation for studying 
caustics is that dark matter, if it annihilates, may do so 
most often in them, since formally the density goes arbi- 
trarily high there. Thus they may greatly enhance observ- 
able signals fro m dark matter (e.g., |Hogan|200"T| [Natarajan] 



fc Sikivie||2008| [Vogelsberger fc White||20li[ ). The behavior 
of caustics in Eulerian space has been extensively studied. 
For example, catastrophe theory allows a rigorous descrip- 
tion and classification of the kinds of caustics and structures 
that can develop in Eulerian space when the manifold folds 
and crosses itself (e.g., [Arnold et al.||1982| |Arnold|2001[ ). 

A Lagrangian view of the cosmic web is difficult to ap- 
ply directly to observed data; however, perhaps methods 
that construct initial conditions from observed data could 



2006 



in principle allow an avenue to do so (e.g. Mohayaee et al, 
Lavaux et al.|20'08 |. And, for example, tracking caus- 



tics in simulations is relevant to the caustic method of esti- 
mating galaxy-cluster masses, which uses the outer caustics 
of cluster haloes as observed in redshift space ( [Diaferio fc 
|Geller|p97| . 

This paper is organized as follows. In Section[2] we very 
briefly review some properties of origami mathematics that 
are applicable in cosmology. In Section |3] we apply these 
properties to cosmology, giving Lagrangian definitions of the 
concepts of 'caustics' and 'streams,' and show some relevant 
measurements from A^'-body simulations. 



2 ORIGAMI MATHEMATICS 

Considering how ancient the art of origami is, the mathemat- 
ics of it have developed relatively recently. A relatively early 
result (dating to 1893) is that paper-folding is a more pow- 
erful mechanism for solving geometric and algebraic prob- 
lems than the classical ruler and compass. For example, in 
principle it can be used to trisect an angle, and to solve 
some types of cubic equations (I Row' 1966' 'Martin"1998' ^HuH 



2006). Then, in the last few decades, significant advances 
were made in algorithmic origami design (e.g., Lang|2003 1. 

Here we focus on some results in 'flat origami' that are 
particularly relevant to large-scale structure. In flat origami, 
folding of a two-dimensional sheet is allowed in three dimen- 
sions, but the result is restricted to lie flat in a plane, i.e. 
it could be squashed between pages in a book without ac- 
quiring any new creases. The class of flat-foldable origami 
is quite large, for example encompassing the famous paper 
crane, similar to the model shown in Fig. [l] 

There are several theorems that have been proven about 
flat origami (e.g., |Hull||1994[ |2002| ). The main flat-origami 
result that we exploit in the present paper is the two- 
colorability of polygons outlined by origami crease lines. 
That is, two colors suffice to color them so that no adja- 
cent polygons share the same color. Here 'adjacent' means 
sharing a crease line (not just a vertex). 

To see why two colors suffice, consider the paper crane 
shown in Fig. [l] Both sides of it are shown, along with its 
appearance when unfolded. Each polygon is colored white or 




Figure 1. Two-coloring of the polygons outlined by creases in an 
origami 'traditional Japanese flapping bird' (similar to a crane). 
Polygons facing 'up' (out of the page in the leftmost diagram, 
with the head facing to the left) are colored white, while polygons 
facing 'down' (into the page in the leftmost diagram) are painted 
gray. The crane has been unfolded for the rightmost diagram. 
Creases are shown as lines here; the polygons outlined by them 
are successfully two-colored. Figure from |HuII| ( |2006^ , courtesy 
Tom Hull. 



gray according to whether the polygon is facing 'up,' i.e. with 
the same orientation as it did initially, or 'down,' if it has 
been flipped over. This uniquely colors each polygon, and 
each crease does indeed divide 'up' from 'down' polygons. 

A work of flat origami can be thought of as a function 
(a continuous piecewise isometry) mapping the unit square 
(the unfolded sheet at right in Fig. [TJ into the plane. Each 
crease produces a reflection, reversing the direction of the 
vector on the paper perpendicular to the crease. The func- 
tion is defined on each polygon by a sequence of these refiec- 
tions. The color in each polygon corresponds to its parity, 
i.e. depending on whether the number of refiections used to 
define the function on that polygon is odd or even. It can 
also be measured locally with the determinant of the matrix 
defining the function on the polygon; we will also use this 
latter definition in the cosmological case below. 

Besides two-colorability, there are other properties that 
fiat-foldable crease patterns have. For example, Maekawa's 
theorem states that in a flat-foldable crease pattern, the 
numbers of 'mountain' and 'valley' creases around a vertex 
differ by two. (A mountain crease becomes folded to form 
an upward-pointing ridge; a valley crease is folded in the op- 
posite way.) There is not an obvious (to the author) exten- 
sion of this result to cosmology, but one may exist. Another 
result for fiat paper origami, Maekawa's theorem, implies 



two-colorability through graph-theory considerations (Hull 



19941. A more difficult property to test is that an arbitrary 



crease-pattern is physically flat-foldable without the paper 



intersecting any folds; this is an NP-complete problem (Bern 
fc Hayes|1996[ ). There are further results that, for instance, 
describe the angles around vertices, but they depend on the 
non-stretchability of the origami sheet, making them inap- 
plicable to the cosmological case. 

Fig. [2] shows a work of fiat origami that bears some 
resemblance to a network of filaments and clusters in cos- 
mology. Its 'voids' are Voronoi cells generated from the 
black points. Voronoi models of large-scale structure are 
good heuristic models of cosmological structure formation 
(e.g., |Icke fc van de Weygaert|p87: Kofman et al.||1990|. 



The present figure corresponds most closely to a Zel'dovich- 
approximation evolution of particle displacements, in which 
structures fold up when expanding voids collide, but over- 
shoot and do not undergo realistic further collapse. Many 
origami tessellations are described, with instructions, by 



Gjerde (20081. 
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berger & White (20111 use this parity measure to explore 



the rococo fine-scale structure of caustics in Eulerian space. 

J(q) gives the volume of a mass element that gets de- 
formed by the tensor, times the parity. The elements of the 
deformation tensor are the components of initial-conditions 
basis vectors transformed to the final conditions, 



J(q) = det 5; 



cm 

dqj 



(1) 



Here is the displacement field. For a mass element with 
Lagrangian position q and Eulerian position x, ^ is defined 
by x(q) = q+^(q). As usual in cosmology, we use comoving 
coordinates. 




Figure 2. A Voronoi origami tessellation that bears some re- 
semblance to cosmological voids, filaments and haloes. If the 
work were unfolded, the polygons outlined by creases would be 
two-colorable according to which way the polygon is facing. E.g. 
'voids' and the topmost polygons in 'haloes' could be colored red, 
and the paper turned upside down within the 'filaments' could 
be colored green. The simple structure of 'filaments' and 'haloes' 
is as in the Zel'dovich approximation; the true universe has more 
complicated folding patterns, shown in subsequent figures. Figure 
by Eric Gjerde ||http: //www. origamitessellat ions . com/\ , used 
with permission. 



3 APPLICATION TO COSMOLOGICAL 
STRUCTURES 

Moving from paper origami to cosmological structure for- 
mation introduces a few changes. The manifold (sheet) has 
three instead of two dimensions. It folds in six dimensions 
(three position, and three velocity) instead of three. The 
sheet also stretches inhomogeneously. Moreover, the metric 
in six dimensions is not necessarily well-defined, although 
one might impose a global metric by e.g. converting veloci- 
ties to positions with the Hubble constant. 

In spite of these difficulties, the local parity of a patch 
in Lagrangian space (the initial conditions) is still easily de- 
fined, in the same way as in fiat origami. Define a first-order 
approximation of the mapping of mass elements (particles) 
from initial to final coordinates as a linear transformation 
(called the deformation tensor) plus a translation. If the spa- 
tial resolution in Lagrangian space (i.e. the mass resolution 
in a simulation) is sufficiently large compared to the scale on 
which Lagrangian space is folding, this linear transformation 
is meaningful. However, deep within a halo, the spatial re- 
lationships between Lagrangian neighbors in a practical A'^- 
body simulation may become essentially random, and this 
condition may not be satisfied. 

We may use the determinant J(q) of this deformation 



tensor to measure the parity (e.g., White & Vogelsberger 



2009); as mentioned above, this definition can be used in 



3.1 Streams and caustics in Lagrangian space 

Caustics and streams are usually considered in Eulerian 
space. A caustic is a fold in projected phase space where 
the density formally goes infinite if particle discreteness is ig- 
nored, and the density is smoothed on arbitrarily fine scales. 
As for streams, at a given Eulerian-space location x, there 
may be many of them; each stream corresponds to a point 
in Lagrangian space that has ended up at x. 

In Lagrangian space, by analogy to an unfolded origami 
sheet, we define a stream as a contiguous three-dimensional 
region with the same parity. We define a caustic as a two- 
dimensional surface separating streams from each other. De- 
fined this way, a caustic indeed corresponds to a fold, since 
the parity swaps if one moves across it. 

By definition, then, space is tessellated by streams that 
are outlined by caustics. The streams are also two-colorable, 
since the parity may take only two values. 

Two-colorability may seem hopelessly academic, and in- 
deed it does not have obvious observational consequences. 
But in fact it greatly restricts the arrangement of streams 
in three-dimensional Lagrangian space. A tessellation in 
greater than two dimensions has in principle no bound on 
its chromatic number (the number of colors required). For 
instance, in relation to the famous four-color theorem (four 
colors suffice to color any planar arrangement of regions), 
Guthrie ([TSSO^ discussed the difficulty of such a result in 
three dimensions. He imagined a set of arbitrarily many 
sticks, each of which touches all others. It is not hard to 
imagine such an arrangement, especially if the sticks are al- 
lowed to fiex. In this case, the chromatic number is bounded 
only by the number of sticks. 

In graph theory, a two-colorable graph is called bipar- 
tite (e.g., 



Chartrand & Zhang 2009 1 . The vertices of our 



bipartite graph are the three-dimensional stream regions, 
and the edges (linking vertices) are the caustic surfaces be- 
tween them. A bipartite graph has many properties; in fact, 



two-dimensional fiat origami, as well. For example, |Vogels-| 



at least one whole book is devoted to the subject (Asratian 
|et al.|[T998[ ). For instance, there is no path that starts and 
ends at the same stream region that consists of an odd num- 
ber of steps from stream to stream. Considering the dual 
graph, in which streams and caustics swap roles, Konig's 
theorem states that the minimum number of colors neces- 
sary for a successful coloring of caustic surfaces equals the 
maximum number of caustic surfaces around a single stream 
region. 

Before continuing on to measurements from simula- 
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tions, we examine our definitions of streams and caustics 
a bit further, and problems that may arise from them. 

A physical problem with our parity definition arises 
from a type of folding that is, in principle, possible in cos- 
mology, but not in paper origami because the paper can- 
not stretch. Cosmological caustics may form in spherical or 
cylindrical collapse, not just planar collapse. Like planar col- 
lapse, spherical collapse reverses parity, but cylindrical col- 
lapse does not; it simply produces a 180° rotation in the two 
axes perpendicular to the cylinder. 

However, in a physically realistic situation, the proba- 
bility that more than one axis will collapse exactly simul- 
taneously is zero. In the Zel'dovich approximation, for ex- 
ample, the deformation tensor will never have two exactly 
equal eigenvalues. 

In a practical A'^-body simulation, the finite temporal 
and mass resolution will enable two- (and three-) axis cross- 
ings to happen between timesteps. In applying our parity 
definition to a simulation, we unfortunately miss some such 
unresolved caustics. 

Another difficulty springs from numerical errors and er- 
rors (e.g. two-body effects) due to particle discreteness in a 
practical simulation. After many dynamical times suffering 
deep within a halo, Lagrangian mass elements may be hope- 
lessly tangled for these reasons, and not just because of the 
doubtlessly plentiful caustics that have formed. Indeed, as 
shown in the figures below, halo regions have fairly random 
parities, and it is not clear whether these parities arise from 
physical caustics, or this numerical 'noise.' 

It may be difficult to visualize how a newly forming 
caustic forming in an already high-density region (already 
with many overlapping streams) behaves in Lagrangian 
space. In Lagrangian space, the caustic surface slices the 
various streams at various angles depending on how each 
stream has been rotated, stretched and reflected. But by 
definition, the new streams are still three-dimensional re- 
gions, divided by two-dimensional caustics. 

3.2 Simulation measurements 

Fig. [3] shows the folding up of a cosmological sheet of 256^ 
particles from a 256''-particle gravitational simulation run to 
the present epoch. Pixels (in Lagrangian coordinates) and 
particles (in Eulerian coordinates) are colored according to 
their ORIGAMI morphologies (left) and parities (right). These 
colorings visually correspond quite well to where structures 
are after folding. There are a few regions that would look 
voidy if they weren't colored, but are colored nonetheless 
(blue at left, and orange at right). These particles have 
crossed others in the direction perpendicular to the page. 

The upper-right panel shows a measurement from an 
A^-body simulation of the quantity we use to determine 
the local parity, J(q) from Eq. [l] This gives the volume 
of the mass element (inversely proportional to its den- 
sity) times its parity. The parity can be read off from 
the color scale. Regions with right-handed parity (as in 
initial conditions) are black or blue; left-handed regions 
are white or orange. Each pixel represents a particle. 
In origami terms, this is the sheet before folding. Note 
that the magnitude is quite small in the cores of halo 
regions, because mass elements shrink considerably in 
high-density halo regions. To estimate the tensor at each 



particle, we use the final-conditions separations between 
the six particles that initially surround the particle, i.e. 

where the indices refer to positions on the Lagrangian 
particle grid. Note that the resolution of this measurement 
is twice the smallest available Lagrangian length scale (the 
interparticle spacing) . 

The simulation sho wn was run to the present epoch us- 
ing the GADGET-2 code ( Springel 2005 I with 256^ particles in 



a box of size cissuming a concordance ACDM set 

of cosmological parameters (fit — 0.04, flm = 0.3, Qa = 0.7, 
h = 0.73, Us — 0.93, as = 0.81). Before running the simu- 
lation, its initial-conditions density field was smoothed with 
a Gaussian of width a = l/i~^Mpc, inhibiting small-scale 
structure formation and clarifying the parity measurement. 
This smoothing reduced as (the dispersion in the overden- 
sity within spheres of radius 8 Mpc) by 0.4 per cent. This 
small-scale power attenuation would happen in the real uni- 
verse if the dark matter had substantial warmth; however, 
warm dark matter would also cause the three-dimensional 
'sheet' it occupies to be smeared out by thermal motions. 

There is some agreement between outer caustics identi- 
fied by ORIGAMI morphology (the boundaries between black 
and non-black regions) and as defined by parity (the outer- 
most boundaries between black/blue and white/orange re- 
gions), but the agreement is not perfect. There are a couple 
of reasons for this. First, the parity is defined on a larger 
grid spacing on the Lagrangian lattice (two interparticle 
spacings, not one). Thus, if only one particle crossed an- 
other, this would be detected by the ORIGAMI criterion but 
not necessarily by the parity criterion. Second, the ORIGAMI 
algorithm detects crossings projecting along the Eulerian 
{x, y, z) axes; in contrast, the parity measurement uses only 
the local vectors between Lagrangian neighbors. So, unlike 
the parity measurement, the ORIGAMI morphology is in prin- 
ciple sensitive to rotations (fortunately, small on cosmologi- 
cal scales). Still, the ORIGAMI halo, filament and wall-finder 
shines here, where the smoothed initial conditions produce 
visually obvious morphologies. 

Fig. [4] is the same figure as Fig. [3] except that the sim- 
ulation was run without smoothing the initial conditions. 
The small-scale structure in Fig. |4] makes the parity map 
(the upper-right panel) much more cluttered from additional 
structure. There are many visible extended streams (patches 
of identical parity), especially voids, but much of the plot 
looks essentially random. Halo particles (comprising most 
of the particles) have likely undergone many shell crossings. 
We again caution, however, that some of this apparent ran- 
domness could be from numerical 'noise.' 

Fig.|5]shows Lagrangian J(q) plots, adding red contours 
at the Zel'dovich-approximation ( |Zerdovich||1970] ) predic- 
tion for the position of the outermost shell crossing. In the 
Zel'dovich approximation, the evolution with time r of the 
volume V (in units of the mean) of a mass element at La- 
grangian coordinate q is given by 



V{ci,r) = l[[l-Kici)Dir)], 



(2) 



a product over the three eigenvalues Xi of the deformation 
tensor, using the linear growth factor D{t). 

For the Zel'dovich prediction in Fig. [S] we computed 
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Figure 3. The folding of a cosmological sheet (top, unfolded; bottom, folded). Quantities were measured from a 256'^ sheet of particles 
from a 256'^-particle ACDM TV-body simulation; The 256^ particles share the same 2-coordinate in the initial-conditions lattice, where 
2 points out of the page. Before running the simulation, the initial conditions were smoothed with a lfe~^Mpc Gaussian window, 
to inhibit small-scale structure formation. Top panels use Lagrangian coordinates, in which each particle is a square pixel in a 256^- 
pixel image. In the bottom panels, particles are shown in their actual present-epoch Eulerian {x, y) coordinates, projecting out the 
z coordinate (in which the slice does have some extent). In left-hand panels, void, wall, filament and halo ORIGAMI morphologies are 
shown in black, blue, yellow and red, respectively. In right-hand panels, particles are colored according to J, i.e. the volume of their 
fluid element times its parity. Black/blue particles have right-handed parity (as in the initial conditions), and white/orange particles 
have swapped, left-handed parity. 



the largest eigenvalue of the deformation tensor (estimated 
as above) at each particle in the initial conditions at redshift 
z = 75. The red curve is the zero-contour of [1 — Amax£'(z ~ 
0)], evolved forward to the present epoch with the ratio of 
growth factors D{z ^ 0)/D{z = 75). 



The agreement is rather good between the Zel'dovich 
prediction and actual outer-caustic locations for the 
smoothed simulation (left). For the unsmoothed simulation, 
the agreement is worse, but still not disasterous. It is not 
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Figure 4. Same as Fig. |3] except measured from a simulation with full initial power, i.e. unsmoothed initial conditions. 



surprising that the Zel'dovich approximation works best at 
finding large-scale caustics. 



4 CONCLUSION 

In this paper, we further elucidate the analogy between 
structure formation in cosmology and origami. We illus- 
trate an insight from origami mathematics that appears to 
be applicable to cosmology: the two-colorability of three- 
dimensional regions (streams) outlined by two-dimensional 
caustic surfaces in the initial conditions. That is, two colors 



suffice to color them such that adjacent regions do not have 
the same color, in the same way that four colors suffice to 
color any planar map. In fact, this result is rather trivial, 
if streams and caustics in Lagrangian space are defined as 
we do. However, a general arrangement of three-dimensional 
regions has no bound on the required number of colors, so 
this is a significant restriction on how dark matter can as- 
semble itself into structures. Two-colorable graphs such as 
this have many properties that may prove useful in under- 
standing structure formation. 

While much is known about the behavior and morphol- 
ogy of caustics and streams in Eulerian space, their behav- 
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Zel'dovich outer-caustic prediction 




Figure 5. Comparison of a Zel'dovich prediction of where outer caustics should form (red contours) to where they actually do 
(boundaries between black and white regions), in the 64-/i~^Mpc simulations. In the left panel, the initial conditions have been 
smoothed, but not in the right panel. Pixels (corresponding to particles) are colored according to the present-epoch deformation tensor 
determinant J(q). The quantities plotted under the red contours are same as in the upper-right panels of Figs, [s] and |4] However, to 
emphasize outer caustics, here we use a different color scale that does not particularly distinguish small negative from small positive 
J. The red contours show the Zel'dovich prediction, i.e. the zero-crossing of the linearly extrapolated smallest eigenvalue of the initial 
deformation tensor. 



ior in Lagrangian space is also of interest. To our knowledge, 
this paper contains the first explicit illustration of the shapes 
of caustics and the streams in Lagrangian space. Many ques- 
tions can be asked about them. How often do they meet in 
one-dimensional curves, and what do these one-dimensional 
curves represent? How often are caustic surfaces isolated? 

Admittedly, two-colorability in itself is lacking in obvi- 
ous observational consequences, but it is closely tied to the 
process of caustic formation, which is of great interest for 
prospective dark-matter direct and indirect detection. As 
cosmological simulations push to smaller and smaller scales, 
departing farther and farther from comfortable linear-regime 
physics, it is useful to know as much mathematically as pos- 
sible (including the present result) about the structures that 
develop. It would be interesting to explore whether enforc- 
ing the properties of caustics and streams found here, or of 
the three-dimensional sheet structure itself, would be useful 
to 'clean' numerical noise away from simulations. The anal- 
ogy to origami is likely to have pedagogical value, as well, 
enabling hands-on paper origami constructions of a cosmic 
web as in Fig. [2] 
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